#include <thrust/detail/device/dereference.h>
#include <cusp/qr.h> /* need rotmat() */

namespace cusp {
namespace device {
namespace block {


/** \brief QR decomposition of a dense matrix.

	Perform qr decomposition of a matrix using one block of threads.

	This implementation uses Givens rotations.

	\note QValues must have been filled with zeros.

	\param values values of matrix Ahat with size num_rows x num_cols
	\param num_rows number of rows in Ahat
	\param num_cols number of columns in Ahat
	\param [out] QValues values of matrix Q with size num_rows x num_rows
	\param [out] RValues values of matrix R with size num_rows x num_cols (same as Ahat)
*/
template<typename MatrixIterator>
__device__
void qr(MatrixIterator values, size_t num_rows, size_t num_cols, MatrixIterator QValues)
{
	using thrust::detail::device::dereference;
	typedef typename MatrixIterator::value_type ValueType;

	__shared__ ValueType tv0[512];
	__shared__ ValueType tv1[512];

	/* make Q a diagonal matrix */
	for (size_t i = threadIdx.x; i < num_rows; i += blockDim.x)
		dereference(QValues + (i * num_rows + i)) = ValueType(1);

	//__syncthreads();

	for (size_t j = 0; j < num_cols; ++j)
	{
		for (size_t i = num_rows - 1; i > j; --i)
		{
			__syncthreads();

			ValueType c, s;
			cusp::detail::rotmat(
				dereference(values + ((i-1) * num_cols + j)),
				dereference(values + ( i    * num_cols + j)),
				c, s);

			__syncthreads();

			for (size_t k = threadIdx.x; k < num_cols; k += blockDim.x)
			{
				tv0[threadIdx.x] = dereference(values + ((i-1) * num_cols + k));
				tv1[threadIdx.x] = dereference(values + ( i    * num_cols + k));

				dereference(values + ((i-1) * num_cols + k)) =  c * tv0[threadIdx.x] + s * tv1[threadIdx.x];
				dereference(values + ( i    * num_cols + k)) = -s * tv0[threadIdx.x] + c * tv1[threadIdx.x];
			}

			for (size_t k = threadIdx.x; k < num_rows; k += blockDim.x)
			{
				tv0[threadIdx.x] = dereference(QValues + ((i-1) * num_rows + k));
				tv1[threadIdx.x] = dereference(QValues + ( i    * num_rows + k));

				dereference(QValues + ((i-1) * num_rows + k)) =  c * tv0[threadIdx.x] + s * tv1[threadIdx.x];
				dereference(QValues + ( i    * num_rows + k)) = -s * tv0[threadIdx.x] + c * tv1[threadIdx.x];
			}
		}
	}
}


} // end namespace block
} // end namespace device
} // end namespace cusp

